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Abstract 

This  project  develops  a  theoretical  model  for  gases  of  bosonic  atoms  at  ultracold,  but 
finite  temperatures.  Under  these  circumstances,  the  gas  can  undergo  a  phase  transition  to  a 
purely  quantum  mechanical  state,  a  Bose-Einstein  condensate  (BEC),  where  the  atoms  cease 
to  behave  like  distinguishable  entities,  and  instead  form  a  single  macroscopic  matter  wave. 
At  exactly  zero  temperature,  all  of  the  atoms  occupy  the  BEC;  at  finite  temperatures, 
a  significant  fraction  of  the  atoms  leave  the  BEC  and  form  a  thermal  cloud.  Thus,  the 
state  of  a  low,  but  finite-temperature  gas  of  bosonic  atoms  involves  the  coexistence  of  a 
BEC  and  a  thermal  cloud.  Further,  the  atoms  can  interact  in  a  variety  of  different  ways, 
which  have  important  consequences  for  the  state  of  the  gas.  We  consider  both  short  range 
contact  interactions  and  dipolar  interactions,  where  the  atoms  interact  via  the  long-range, 
anisotropic  dipole-dipole  force.  We  develop  this  model  in  both  three-  and  two-dimensional 
geometries. 

When  the  dipolar  gas  is  confined  to  a  nearly  two-dimensional  geometry,  the  BEC  can 
exhibit  an  unusual  excitation  spectrum,  where  certain  excitations  that  would  usually  be 
highly  energetic  are  instead  only  weakly  energetic.  These  excitations  are  called  rotons. 
Because  it  costs  minimal  energy  for  an  atom  to  leave  the  BEC  and  be  promoted  to  a 
roton  excitation,  modest  temperatures  can  result  in  fairly  large  thermal  cloud  occupations. 
This  project  describes  the  relative  occupations  of  the  BEC  and  the  thermal  cloud  at 
given  temperatures,  and  when  roton  excitations  are  present.  Further,  BECs  in  nearly 
2D  geometries  take  the  form  of  quasi-condensates,  or  BECs  with  finite  spatial  extent. 
Quasi-condensates  behave  like  BECs  on  shorter  length  scales,  but  not  on  longer  length 
scales.  The  project  incorporates  the  presence  of  a  quasi-condensate  into  the  theoretical 
model,  to  more  accurately  describe  the  state  of  the  system  at  hand. 


Keywords:  Bose  Einstein  condensation,  ultra  cold  physics,  condensed  matter,  dipoles 
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I.  BACKGROUND 

Quantum  mechanics,  the  science  of  tiny  objects,  plays  an  increasingly  significant  role  in 
modern  technologies,  and  is  of  inherent  interest  due  to  its  many  peculiarities.  For  example, 
quantum  mechanics  predicts  that  atoms  possess  wave-like  properties,  such  as  a  wavelength, 
and  can  be  classified  as  either  ’’fermions”  or  ’’bosons.”  Thus,  every  particle  has  a  wavelength, 
yet  also  behaves  like  an  individual  particle  with  either  bosonic  or  fermionic  properties.  Ad¬ 
ditionally,  both  wave  and  particle  treatments  are  necessary  to  fully  understand  the  behavior 
of  quantum  objects.  Concerning  the  particle  classifications,  examples  of  fermions  are  the 
matter-particles  that  everyone  learns  about  in  chemistry  class:  neutrons,  protons,  and  elec¬ 
trons.  Bosons  are  generally  seen  as  force-carrier  particles,  or  particles  that  propagate  the 
four  fundamental  forces  (photons  are  bosons,  for  example).  However,  certain  atoms,  such  as 
Helium-4  (4He),  which  is  composed  entirely  of  fermions,  acts  like  a  bosons.  It  turns  out  that 
this  is  the  case  whenever  an  object  like  an  atom  is  composed  of  an  even  number  of  fermionic 
constituents.  These  atoms  are  called  “boson-like,”  [1],  though  here  we  will  simply  refer  to 
them  as  bosons.  Whereas  fermions  are  not  allowed  to  occupy  the  same  point  in  space  due  to 
the  Pauli  Exclusion  Principle,  bosons  can  occupy  the  same  point  in  space.  This  is  the  major 
distinction  between  them.  Quantum  mechanically  speaking,  this  corresponds  to  fermions 
not  being  able  to  occupy  the  same  quantum  state,  with  bosons  being  able  to  occupy  the 
same  quantum  state.  In  fact,  an  unlimited  number  of  bosons  can  occupy  a  given  quantum 
state.  If  electrons  were  bosons,  they  could  then  occupy  the  same  orbitals  in  atoms,  and  the 
chemistry  of  our  world  would  certainly  be  a  lot  different. 

The  wavelength  of  atoms  is  qualitatively  determined  by  the  De  Broglie  wavelength  A, 
which  is  inversely  proportional  to  the  atom’s  momentum  p,  X  =  h/p ,  where  h  =  1.055  x 
10-34  J  •  s  is  the  reduced  Planck  constant.  In  a  gas,  only  when  the  inter-particle  spacing 
approaches  the  order  of  the  de  Broglie  wavelength  do  quantum  mechanical  effects  begin  to 
matter.  Thus,  at  relatively  high  temperatures  when  the  average  momentum  of  the  atoms  is 
large,  the  properties  of  an  atomic  gas  do  not  depend  on  whether  the  atoms  are  fermions  or 
bosons;  i.e.  the  high  temperature  gas  is  classical.  As  the  temperature  drops  to  the  ultra-cold 
regime,  the  kinetic  energy  each  particle  has  decreases,  and  the  particle  wavelengths  grow 
and  begin  to  overlap.  If  the  atoms  are  bosons,  they  form  a  Bose-Einstein  condensate  (BEC). 
The  atoms  can  no  longer  be  described  individually  and  instead  behave  like  a  giant  matter 
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wave. 

The  Bose-Einstein  condensate  received  its  name  from  two  instrumental  physicists  and 
mathematicians  that  heavily  contributed  to  its  finding.  Satyendra  Nath  Bose  earned  a 
bachelor’s  and  a  master’s  degree  at  the  age  of  twenty  and  then  began  translating  Albert 
Einstein’s  papers  on  general  and  special  relativity  in  order  to  better  understand  relativity. 
As  he  progressed  in  his  studies,  he  began  to  delve  into  quantum  mechanics  more.  Eventually, 
Bose’s  contributions  to  quantum  physics  would  be  so  significant  Paul  Dirac,  one  of  the 
greatest  physicists  of  the  20th  century,  would  name  the  boson  particle  in  his  honor.  In  1924, 
Bose  derived  Planck’s  quantum  radiation  law  with  no  reference  to  classical  physics.  In  order 
to  be  published,  he  sent  the  paper  to  Einstein  to  translate  it. 

The  other  contributing  scientist  is  world  famous  physicist  Albert  Einstein,  maybe  best 
known  for  his  mass-energy  equivalence  equation,  but  he  made  many  more  extensive  con¬ 
tributions  to  physics.  His  works  include,  but  are  not  limited  to  special  relativity,  general 
relativity,  evolving  quantum  theory,  and  the  explanation  of  the  photoelectric  effect.  How¬ 
ever,  his  significant  contribution  to  Bose-Einstein  Condensates  lies  in  response  to  receiving 
Bose’s  work  on  quantum  statistics  of  light  quanta  (photons). 

Einstein  thoroughly  enjoyed  the  paper  so  he  used  his  influence  to  publish  the  paper  and 
personally  translated  it  to  German.  He  then  extended  Bose’s  idea  to  atoms,  which  led  to  the 
creation  of  Bose-Einstein  statistics.  Bosons  are  indistinguishable  from  each  other  and  can 
occupy  the  same  energy  levels.  Therefore,  we  count  indistinguishable  bosons  and  fermions 
much  differently.  Indistinguishable  versus  distinguishable  statistics  can  be  eloquently  shown 
through  coin  flipping.  Say  we  want  to  count  the  probability  of  flipping  two  coins  and  getting 
a  heads  (H)  and  a  tails  (T).  Common  sense  dictates  that  there  are  four  possibilities  (or  mi¬ 
crostates),  HH,  HT,  TH,  and  TT  because  you  can  individually  identify  which  coin  is  which. 
Maybe  the  two  coins  are  a  nickel  and  a  quarter  or  quarters  of  different  years.  The  important 
part  is  that  you  can  distinguish  between  them  and  therefore  can  distinguish  separate  mi¬ 
crostates.  With  these  distinguishable  coins,  a  tails  and  heads  combination  would  have  two 
out  of  four  total  microstates  or  a  |  chance  of  occurring.  If  the  coins  were  indistinguishable, 
the  statistics  change.  Consider  taking  two  newly  minted  identical  quarters.  Since  there  is  no 
possible  way  to  tell  the  coins  apart,  HT  and  TH  are  the  same  microstate.  Because  of  this, 
the  microstates  would  be  HH,  HT,  and  TT.  HT  would  have  |  chance  of  occurring.  This  dif¬ 
ference  accounts  for  the  seemingly  innocuous,  but  crucial  difference  in  the  statistics  used  to 
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analyze  bosons.  Following  are  the  different  equations  for  Maxwell-Boltzmann,  Fermi-Dirac 
and  Bose-Einstein  statistics 


Maxwell  —  Boltzmann  :  m  =  — — -  '  — 

e{ei-it)/KbT 


(1) 


Fermi  —  Dirac  :  n,  = 


9i 


(2) 


Bose  —  Einstein  :  rq  = 


9i 


(3) 


In  Eq  1,  Eq  2,  and  Eq  3,  rti  describes  the  number  of  particles  in  a  given  state,  i.  The 
state  is  given  by  the  energy  level  associated  with  it  et.  kb  is  the  Boltzmann  constant  and  T 
stands  for  the  temperature  at  which  the  state  is  being  measured,  /i  is  a  Lagrange  multiplier 
used  to  fix  particle  number.  If  we  are  modeling  photons,  we  do  not  require  constant  particle 
numbers  so  we  can  ignore  /i.  Finally,  r/,  is  the  degeneracy  associated  with  the  selected  state. 
If  two  states,  e2  and  e3  share  the  same  energy  level,  e2,  the  states  are  called  degenerate  and 
gi  will  account  for  it. 

Equation  1,  Maxwell-Boltzmann  statistics,  describes  the  distribution  of  non-interacting 
particles  in  thermal  equilibrium  at  temperatures  high  enough  (or  density  low  enough)  to 
ignore  quantum  effects.  We  use  these  statistics  to  model  a  gas  (whether  indistinguishable  or 
distinguishable)  at  high  temperatures  where  the  gas  behaves  classically.  Equation  2,  Fermi- 
Dirac  statistics,  describes  the  distribution  of  indistinguishable  fermions,  particles  that  obey 
the  Pauli-Exclusion  Principle.  These  statistics  cover  the  quantum  realm,  but  specifically 
with  regards  to  fermions.  As  previously  stated,  the  inability  of  fermions  to  occupy  the  same 
quantum  space  requires  we  treat  them  differently.  Finally,  equation  3,  Bose  Einstein  statis¬ 
tics,  describes  the  energy  level  occupations  of  indistinguishable  bosons.  [2]  These  statistics 
cover  the  quantum  realm  with  regards  to  bosons. 

In  1924,  Einstein  predicted,  based  purely  on  Eq.  3  and  Bose-Einstein  statistics,  a  unique 
phenomenon  occurs  at  low  temperature.  [3]  As  temperature  decreases,  a  particle  will  lose 
energy,  leaking  it  to  the  environment.  Near  absolute  zero,  this  particle  will  most  likely 
occupy  its  ground  state,  or  least  energetic  state.  When  a  collection  of  bosons  all  occupy  the 


ground  state  simultaneously,  a  new  state  of  matter  called  the  BEC  appears.  In  Einstein’s 
era,  however,  there  was  no  way  to  cool  anything  to  the  temperatures  required  to  observe 
physical  BECs,  so  physical  exploration  was  tabled  until  adequate  cooling  techniques  existed. 

As  time  progressed,  more  research  could  be  made  towards  the  ultimate  goal  of  creating 
this  new  state  of  matter.  Superfluidity  was  first  observed  in  liquid  Helium-4  in  1937,  at 
temperatures  below  T=2.17  Kelvin.  In  1938,  Fritz  London  proposed  BECs  as  a  mecha¬ 
nism  for  superfluidity.  Superfluidity,  possibly  the  most  interesting  trait  of  BECs,  is  when  a 
liquid/current  flows  with  absolutely  zero  viscosity  below  a  certain  critical  temperature.  [4] 
Interestingly,  this  is  the  underlying  mechanism  for  superconductivity  in  many  metal  and 
ceramic  materials.  If  two  electrons  pair  together,  they  act  like  a  boson.  Once  you  have  a 
collection  of  bosons  (electron  pairs),  they  can  condense  and  achieve  superfluid  properties. 
However,  you  need  multiple  bosons  to  achieve  a  BEC.  A  single  boson  cannot  become  a  BEC, 
but  many  pairs  of  electrons  will  be  able  to  condense.  It  is  believed  that  the  superfluidity 
in  liquid  Helium-4  is  due  to  the  presence  of  a  BEC  in  the  liquid,  though  this  is  difficult  to 
show  directly.  However,  not  all  BECs  produce  superfluids,  and  not  all  superfluidity  stems 
from  BECs.  Superfluidity  remains  desirable,  because  it  would  eliminate  energy  losses  in 
fluid  systems  and  electrical  systems. 

Despite  continued  predictions,  Bose-Einstein  Condensates  remained  elusive  due  to  the 
technological  restraint  of  lowering  temperatures  down  to  near  absolute  zero.  In  1978,  laser 
cooling  was  invented  and  the  ability  to  reach  previously  unfathomable  temperatures,  within 
nanoKelvin  of  absolute  zero,  became  realistic.  As  the  technology  developed,  this  type  of 
cooling  became  readily  available  and  practical.  Laser  cooling  is  successful  because  photons 
have  momentum.  If  a  photon  is  traveling  in  the  opposite  direction  of  an  atom  and  is 
absorbed,  the  atom  will  slow  (Imagine  a  collision  between  a  bowling  ball  rolling  down  an 
alley  and  a  ping  pong  ball  shot  in  the  opposite  direction.  The  bowling  ball  slows,  but 
not  significantly  enough  to  be  noticed).  Then,  the  atom  emits  the  photon  in  a  random 
direction.  On  average,  this  will  additionally  decrease  the  momentum  of  the  atom.  This 
can  be  imagined  as  the  bowling  ball  absorbed  the  ping  pong  ball  and  randomly  ejecting 
it.  In  half  of  the  possible  ejections,  the  bowling  ball  will  slow  even  more,  albeit  minimally. 
However,  after  enough  absorptions  (nearly  20,000  for  a  sodium  atom),  the  atom  will  near 
absolute  zero.  [5]  Because  a  laser  can  induce  up  to  107  absorptions  per  second,  this  process 
can  occur  in  milliseconds. 


9 


With  this  new  technology,  scientists  could  finally  examine  Einstein’s  predictions  about  the 
BEC.  On  June  5,  1995,  University  of  Colorado  Boulder  became  the  first  to  synthesize  a  Bose- 
Einstein  Condensate  with  Rubidium  (Rb)  atoms  at  near  170  nK.  [6]  That’s  1.7  x  10-7  Kelvin 
above  absolute  zero,  which  was  the  coldest  known  temperature  in  the  universe!  Within  a 
few  months,  MIT  researchers  also  synthesized  a  BEC,  but  with  Sodium  (Na)  atoms.  [7]  Both 
teams  received  the  2001  Nobel  Prize  for  their  foray  into  a  new  section  of  physics.  [3] 

The  study  of  BECs  recently  underwent  even  more  dynamic  and  interesting  developments. 
In  1998,  just  three  years  after  a  BEC  was  first  realized,  Lene  Hau  of  Harvard  managed  to 
slow  light  down  to  17  meters  per  second  by  trapping  it  within  a  BEC.  Within  the  following 
years,  Hau  and  her  team  managed  to  not  only  slow  light  down,  but  capture  it  outright. 
They  stopped  a  pulse  of  light  within  a  BEC,  released  it,  and  then  recaptured  the  light 
packet  within  another  BEC.  This  raises  the  possibility  of  new  developments  in  light  based 
communications,  optical  data  storage,  and  quantum  computing. 

Today,  BECs  can  be  made  with  a  variety  of  atomic  species,  including  Chromium  (Cr)  [8] 
and  the  rare-earth  (lanthanide)  atoms  Erbium  (Er)  [9]  and  Dysprosium  (Dy)  [10],  which 
possess  large  magnetic  dipole  moments.  Further,  BECs  are  used  for  a  variety  of  purposes, 
including  the  investigation  of  matter’s  fundamental  low-temperature  properties.  For  ex¬ 
ample,  superfluid  flow  and  quantized  circulation  (quantum  vortices)  were  directly  observed 
in  experiments  with  BECs.  [11]  [12]  The  National  Institute  of  Standards  and  Technology 
in  Gaithersburg,  Maryland  carried  out  many  of  these  experiments.  However,  the  experi¬ 
ments  involved  atoms  that  possess  negligible  dipole  moments  (like  the  alkali  atoms).  Today, 
there  is  a  strong  effort  to  learn  more  about  the  nature  of  BECs  of  dipolar  atoms,  the  so- 
called  dipolar  BECs.  Unlike  most  atoms,  which  interact  via  the  short-range  Van  der  Waals 
force,  dipolar  atoms  interact  via  the  dipole-dipole  force,  which  is  anisotropic  and  long-range. 
When  oriented  side-by-side,  two  dipoles  repel  each  other,  and  when  oriented  head-to-tail, 
they  attract.  For  this  reason,  experiments  with  dipolar  BECs  usually  work  in  nearly  two  di¬ 
mensional  (2D)  geometries,  with  the  dipoles  polarized  perpendicular  to  the  2D  plane.  This 
results  in  a  configuration  where  the  dipoles  are  oriented  side-by-side,  and  exert  repulsive 
forces  on  each  other,  which  helps  to  stabilize  the  BEC.  The  presence  of  the  dipole-dipole 
force  has  been  shown  to  result  in  interesting  phenomena  in  dipolar  BECs,  and  a  growing 
body  of  theoretical  work  provides  motivation  for  future  experiments.  [13] 

The  attractive  interaction  results  in  a  unique  excitation  called  a  roton  excitation.  A  roton 
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excitation  allows  for  a  particle  with  high  momentum  and  short  wavelength  to  be  excited  out 
of  the  condensate  for  very  low  energy.  This  is  atypical.  In  order  to  explore  how  rotons  affect 
a  BEC  and  its  superfluid  properties,  we  first  need  to  understand  the  mechanisms  already 
in  place  to  describe  BECs.  From  a  non-interacting  case,  to  contact  interactions,  and  finally 
to  dipolar  interactions,  the  following  section  will  provide  the  knowledge  to  understand  the 
functioning  of  the  BEC. 


II.  THEORY 

A.  Non-Interacting  Bose  Gas 

The  original  phenomenon  of  Bose  Einstein  condensation  focused  on  a  non-interacting  gas. 
Einstein  predicted  Bose  Einstein  condensation  because  it  is  a  purely  statistical  phenomenon. 
This  means,  that  a  zero  interaction  BEC  exists  only  due  to  quantum  statistics  and  not  to 
interactions  between  molecules,  intra  molecular  interactions,  or  the  nature  of  the  particle.  [4] 
For  a  BEC  to  exist,  there  must  be  a  multitude  of  particles.  A  single  particle  cannot  become 
a  BEC.  Being  a  purely  statistical  phenomenon  requires  large  numbers  of  particles,  despite 
that  each  individual  particle  has  no  idea  any  other  particles  are  around  it.  As  mentioned  in 
equation  3,  the  Bose  Distribution  governs  the  statistical  occupation  of  the  condensate, 


Hi  e(ei—/j,)/kbT  _  x  W 

Here,  denotes  the  energy  of  a  single  particle  state,  p  is  a  chemical  potential  that  acts 
as  a  Lagrange  multiplier  to  fix  particle  number  in  the  Grand  Canonical  Ensemble.  In  order 
to  achieve  this,  the  condition  is  applied  that  the  sum  of  all  individual  occupancies  must  sum 
to  the  total  particle  number,  /p,  is  the  Boltzmann  constant  and  T  is  temperature.  We  fully 
calculate  the  Bose-Einstein  Distribution  because  we  know  the  energy  levels  in  free  space  for 
non-interacting  particles.  When  interactions  appear  in  our  condensate,  calculating  becomes 
significantly  more  difficult,  as  we  will  discuss  later.  Based  solely  on  quantum  statistics,  the 
Bose  Distribution  dictates  the  occupancy  of  each  state,  ground  and  excited.  [4] 

At  high  temperatures,  these  quantum  statistics  are  meaningless.  Any  large  temperature 
precludes  the  ability  for  a  condensate  to  form  and  the  statistics  that  govern  it  become  the 
Maxwell  Boltzmann  statistics  mentioned  in  Eq  1  [2], 
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However,  Bose-Einstein  statistics  are  required  near  the  critical  temperature  because  quan¬ 
tum  statistics  fully  control  the  arrangement  of  bosons  at  such  low  temperatures.  The  critical 
temperature,  defined  as  the  highest  temperature  at  which  the  macroscopic  occupation  of  the 
lowest  energy  state  appears,  is  the  lowest  temperature  in  which  all  atoms  are  in  excited  states 
(no  atoms  in  the  ground  state).  The  maximum  values  of  particles  in  excited  states  results 
in  a  /i  of  zero.  The  equation  for  number  of  particles  in  an  excited  state  is  given  as  followed 
Nex  =  f  dk/(2n)3fk.  /*,  is  the  same  as  the  Bose-Einstein  distribution  mentioned  earlier  in 
Eq  3,  we  simply  replace  rii  with  / *.  (k  being  the  new  momentum  quantum  number).  The 
integral  is  then  rewritten  in  terms  of  x  —  e/kbTcrit  and  evaluated.  After  changing  bases  back 
and  solving,  the  calculations  result  in  [14] 


Tcrit 


2nh2 


N 


mkb  L^C(3/2) 


2/3 


(5) 


If  we  rewrite  the  result  ofEq5,  we  get  Nex  =  N(T /Tcrit)a .  We  find  N0/N  =  [1— (T /Tcrit)a] 
where  a  is  3/2  for  particles  in  a  3  dimensional  (3D)  box  of  volume  V.  [14]  This  result  provides 
for  the  graphing  of  the  condensate  fraction  versus  the  temperature  critical  temperature 
fraction  which  will  result  in  Fig.  1.  The  non-interacting  BEC  is  a  purely  quantum  statistical 
phenomenon  able  to  be  efficiently  and  gracefully  described  using  the  above  equations. 

The  Schrodinger  equation  is  the  only  way  to  describe  a  quantum  object.  This  equation 
— /i2/2m2V2T  +  I/T  =  AT  can  describe  the  motion  and  energy  levels  of  any  particle.  [4] 
Schrodinger’s  equation  can  be  written  as  =  AT,;,  because  the  left  side  of  Schrodinger’s 
equation  is  called  the  Hamiltonian  (H).  h  is  a  constant  and  m  is  the  mass  of  the  particle.  In 
order  to  simplify  future  equations,  we  set  h  =  m  =  1.  T  is  the  wave  equation  that  describes 
the  motion  of  the  particle.  We  can  solve  the  differential  equation  for  various  T’s  depending 
upon  the  potential  energy  function,  U,  used  in  the  equation.  For  a  particle  in  free  space 
(I/=0),  we  find  that  the  type  of  waves  that  satisfy  it  are  plane  waves  that  take  the  form, 
T  =  Aelkx~lU}t.  Although  there  are  two  types  of  Schrodinger’s  equations,  time  independent 
and  time  dependent,  we  use  the  time  independent  equation  because  it  gives  us  energies  which 
we  can  plug  into  the  Bose-Einstein  distribution.  Without  the  time  component,  plane  waves 
take  the  form  T  =  Aelkx  with  A  =  1/LD^2 .  D  is  the  dimensionality  of  our  system.  L  is  the 
length  scale  over  which  the  quantum  object  is  being  described.  If  confined  by  an  infinite 
potential  well,  L  is  the  distance  between  the  two  walls  of  the  well.  In  order  to  find  the 
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FIG.  1:  Condensate  fraction  of  a  non-interacting  BEC.  At  0  Kelvin,  every  atom  is  within  the 
condensate.  As  temperature  increases,  the  condensate  fraction  begins  to  drop,  more  slowly  at  first, 
but  gradually  sharpens  as  it  reaches  the  critical  temperature.  Beyond  the  critical  temperature, 
every  atom  occupies  the  thermal  gas. 

energy  levels  the  particle  can  inhabit,  we  operate  on  the  wave  equation  with  the  quantum 
Hamiltonian  operator.  The  return  of  this  operation  gives  the  eigenvalues,  or  energy  levels,  of 
6/;  =  ti2k2/2m.  k  is  a  quantum  number  describing  the  state  level  that  the  particle  occupies. 
At  low  energy  in  a  box  of  length  LD ,  the  particle  can  only  exist  in  discrete  energy  levels.  [14] 
At  higher  energies  with  no  potential,  the  energy  levels  are  not  defined. 

Since  non-interacting  Bose-Einstein  condensation  is  purely  statistical,  we  require  large 
numbers  of  particles.  For  a  BEC  without  interactions,  each  particle  feels  no  effect  from 
any  other  particle.  Therefore,  we  treat  the  condensate  and  surrounding  thermal  cloud  as  a 
collection  of  individual  particles.  For  a  single  particle,  the  Schroclinger  Equation  gives  the 
energy  which  we  can  plug  into  the  Bose  Distribution  for  a  given  temperature. 
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B.  The  Challenge  Posed  By  Interactions 

If  we  have  several  interacting  atoms  to  describe,  the  Hamiltonian  becomes  more  compli¬ 
cated  H  =  2m^  +  Si  ^(d)  +  \  Sij  ^ ('A  —  h})-  I11  order  to  account  for  the  interactions 

between  particles,  we  add  a  term  to  our  potential  to  describe  particle-particle  interactions. 
An  example  of  what  this  terms  would  look  like  is  V (ft  —  fj)  =  Sr  4^eVl^  —  D |  hi  the 
case  of  electrons  repelling  via  Coulomb  interactions.  For  one  particle,  our  solution  remains 
unchanged.  For  two  particles,  the  solution  becomes  more  difficult  but  still  reasonable.  When 
our  particles  do  not  interact,  our  wave  equation  is  a  product  of  single  particle  wave  func¬ 
tions.  When  our  particles  do  interact,  each  particle  must  remain  symmetric  upon  exchange 
with  any  other  particle.  As  we  increase  particle  number  drastically,  keeping  track  of  each 
particle’s  interactions  with  every  other  one  is  absurdly  difficult. 

This  difficulty  is  compounded  by  the  difference  between  bosons  and  fermions  and  how 
they  cause  the  wave  functions  to  alter.  When  we  have  the  requirement  that  particles  are 
indistinguishable  like  our  particles,  we  must  ensure  that  the  probability  density  of  the  wave 
function  remains  symmetric  upon  exchange.  If  the  probability  density  differed  when  swap¬ 
ping  two  particles,  then  the  particles  are  distinguishable.  There  are  two  methods  of  con¬ 
structing  a  wave  function  to  describe  indistinguishable  particles.  The  first  method  is  by 
ensuring  that  when  swapping  two  particle  positions,  the  sign  of  the  wave  equation  remains 
the  same  d^aq,^)  =  ^(oq, aq).  This  is  called  a  symmetric  wave  function.  The  second 
method  is  by  ensuring  that  when  two  particles  swap  positions,  a  minus  sign  is  introduced 
H/(xi,X2)  =  — T(x2,xi).  This  is  called  an  antisymmetric  wave  function.  When  a  particle  is 
under  symmetric  exchange,  it  is  a  boson.  Likewise,  particles  which  are  antisymmetric  under 
exchange  are  fermions.  Since  we  model  bosons,  all  our  particles  will  be  symmetric  under 
exchange. 

The  Schrodinger  equation  is  capable  of  fully  describing  a  single  particle.  With  two  parti¬ 
cles,  we  treat  the  system  in  its  center  of  mass  frame  (effectively  framing  it  as  a  one  particle 
problem).  As  we  transition  to  a  many  body  system,  particle  number  ranging  from  3  to 
infinity,  the  equations  become  increasingly  difficult  to  solve.  We  can  numerically  solve  many 
body  problems,  but  computations  become  exponentially  more  difficult  the  more  particles  we 
encounter.  Every  particle  has  3  spatial  degrees  of  freedom  (in  3D)  so  a  system  of  N  particles 
has  3N  degrees  of  freedoms.  In  addition,  a  symmetric  wave  function  requires  that  the  wave 
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function  remain  unchanged  through  particle  exchange.  This  means  a  system  of  two  particles 
would  look  like  T  =  ^('0i(^i)V,2(r2)  +  0 i(r2)02(ri))-  A  system  of  3  particles  would  have 
6  terms  in  it.  As  we  increase  the  number  of  particles  in  the  system,  the  number  of  terms 
increases  as  N\.  Even  if  N  is  just  20,  the  number  of  terms  in  our  wave  function  we  need 
to  keep  track  of  is  over  2.4  x  1018.  The  systems  of  particles  we  intend  to  work  with  will  be 
much  more  vast  than  twenty  so  a  different  method  of  computation  is  required.  We  describe 
our  method  for  overcoming  these  difficulties,  ’’second  quantization”,  in  Section  II C. 

1.  Contact  Interactions 

Although  we  will  focus  on  dipolar  interactions,  most  atoms  exhibit  contact  interactions. 
Van  der  Waals  interactions,  the  normal  mode  for  most  atoms,  model  contact  interactions 
well.  Particles  with  contact  interactions  can  be  thought  of  as  pool  balls.  Our  particles 
move  and  if  they  make  contact  with  another  particle,  they  exchange  energy  and  momentum. 
The  particles  then  move  off  in  different  directions,  until  they  interact  with  another  particle. 
Without  contact  interactions,  the  particles  would  pass  right  through  each  other,  completely 
unaware  of  the  other’s  existence. 

We  model  contact  interactions  using  a  delta  function  and  a  contact-interaction  coupling, 
g,  which  we  can  adjust  to  model  stronger  and  weaker  interactions.  Physically,  g  corresponds 
to  a  ratio  of  the  scattering  length  of  the  particles  to  the  inter  particle  spacing.  For  our 
dilute  gas,  g  —  1  models  the  strongest  interaction  strength  we  will  use  while  g  =  0  models 
no  interactions  at  all.  The  potential  itself  manifests  as  V(r  —  r')  =  gd(r  —  r'). 

2.  Dipolar  Interactions 

Magnetic  moments  in  an  atom  come  from  the  total  accumulation  of  electron  spin  added 
with  orbital  angular  momentum  of  that  specific  atom.  The  magnetic  dipole  moment  can 
be  calculated  from  the  total  angular  momentum  of  the  atom  via  the  following  equation 
^Atom  =  fill1  R \/ j (j  +  !)■  Ill  this  equation,  j  is  the  total  angular  momentum  and  gj  and  /jlb 
are  constants. 

For  dipoles  with  dipole  moments  of  magnitude  d  all  polarized  in  the  same  direction,  d  =  z 
describes  a  dipole  moment  polarized  along  the  z  axis.  [15]  The  dipole  dipole  interaction 
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potential  is  given  by 


Vd(r)  =  d 2 


1  —  3  cos  9 


r*/  1 3 


(6) 


6  is  the  angle  between  the  dipole  polarization  vector  d  and  the  vector  separating  the 
dipoles  \  f  —  r'\.  [16]  Vectors  r  and  r'  give  the  position  of  the  two  dipoles  and  d  is  the  dipole 
moment.  These  can  range  from  zero  to  completely  non-interacting  (in  terms  of  dipoles)  to 
upwards  of  10  in  dimensionless  units.  This  constant  is  determined  by  the  physical  dipole 
moments  of  the  atom. 

The  dipole-dipole  interaction  can  be  attractive  or  repulsive.  When  the  dipoles  align  hori¬ 
zontally,  they  repel  each  other.  When  the  dipoles  instead  align  vertically,  they  attract.  This 
causes  problems  for  a  BEC  because  particles  that  are  too  attractive  will  acquire  too  much 
kinetic  energy  exciting  them  out  of  the  BEC.  In  order  to  avoid  this,  we  confine  our  particles 
to  a  2D  geometry  so  they  only  sample  the  repulsive  potential  of  each  other.  However,  as  the 
momentum  of  the  particles  increases,  their  size  nears  the  same  order  as  the  trap.  This  allows 
the  bosons  to  begin  to  sample  the  attractive  potential  of  dipolar  interactions  as  oppose  to 
solely  the  repulsive  component. 


C.  Second  Quantization 

An  alternative  method  to  perform  calculations  without  the  Schrodinger  is  by  using  2nd 
Quantization.  Our  calculations  for  2nd  Quantization  are  shown  in  appendix  A.  2nd  Quanti¬ 
zation  represents  the  quantum  many  body  states  as  part  of  a  Fock  state  basis.  Fock  states 
are  a  form  of  a  Hilbert  space.  Hilbert  space  is  an  infinite  orthonormal  basis  set.  For  ex¬ 
ample,  Cartesian  coordinates  scale  where  each  new  dimension  becomes  orthonormal  to  all 
previous  dimensions.  As  we  gain  more  dimensions,  the  orthogonality  of  our  basis  remains 
the  same.  Fock  states  are  a  quantum  mechanical  states  that  have  well  defined  number  of 
particles  within  them.  For  example  a  state  with  2  particles  in  the  ground  state,  none  in  the 
first  excited  state,  three  in  the  second  excited  state,  and  seven  in  the  third  excited  state  will 
appear  as  xi  —  |2,  0,  3,  7, ...  >.  In  order  to  change  the  number  of  particles  in  a  certain  state, 
we  operate  on  it  withe  a  creation  or  an  annihilation  operator.  dk=2Xi  —  \/n|2,  0, 4,  7, ...  > 
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where  n  is  the  number  of  states  in  the  wave  equation. 

In  order  to  convert  from  our  Hilbert  space  wave  function  to  Fock  space  wave  functions 
(Fock  space  is  a  specific  type  of  Hilbert  space),  we  expand  'k  and  dA  as  a  function  of  an 
orthonormal  basis  with  creation  and  annihilation  operators.  [3] 


^(f)  =  ^y 


(7) 


'k(r)  = 


\/~Ld 


£■ 


iq-r ; 


(8) 


Creation  and  annihilation  operators  do  exactly  what  their  names  suggest.  One  creates  a 
particle  in  a  specific  state  and  one  annihilates  a  particle  in  a  specific  state.  [3]  Creation  and 
annihilation  operators  are  referred  to  as  raising  and  lowering  operators  as  well. 


\Nk=o’  y  ,• 
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We  substitute  these  new  'ks  into  our  non-interacting  Hamiltonian, 


H 


(10) 


and  using  several  commutation  relations  as  seen  in  Appendix  A,  arrive  at  a  new  Fock  space 
Hamiltonian. 


tf1  =  £asaey  <n) 

k 

We  use  ensure  our  Hamiltonian  maintains  symmetry  by  relying  on  the  commutation 
relations.  Our  system  is  made  up  entirely  of  bosons  so  we  require  a  symmetric  wave  function, 
despite  introducing  second  quantization. 

In  addition,  we  replace  k2  / 2  with  eg;  the  single  particle  energy  term,  eg  describes  the 
kinetic  energy  of  a  single  particle.  Once  summed,  we  can  find  the  kinetic  energy  of  our  whole 
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system.  As  previously  mentioned,  /i  is  the  chemical  potential  and  serves  as  a  Lagrange 
multiplier  to  fix  particle  number.  It  appears  here  from  writing  our  Hamiltonian  in  the 
grand-canonical  ensemble, 


H  =  -»)  ■  (12) 

k 

We  perform  the  same  type  of  operation  for  the  interaction  potential  and  the  interaction 
Hamiltonian. 


Hint  —  2 J^D  V.  |  J^k^P  (12) 

k,p,q 

Here,  V (q)  is  the  Fourier  transform  of  V (f)  into  momentum  space.  This  term  describes 
the  exchange  of  momentum  between  two  particles  when  they  interact.  We’ve  begun  with 
modeling  a  general  contact  interaction  potential,  but  eventually  will  focus  on  magnetic 
dipoles  in  a  quasi-2D  case. 


D.  Mean  Field  Approximation 

Despite  the  work  put  into  second  quantization,  our  Hamiltonian  is  still  not  in  a  practical 
form.  We  need  the  Hamiltonian  to  be  diagonal  and  quadratic  in  order  to  find  its  eigenvalues 
(energy  of  the  system).  Our  first  step  in  simplifying  our  equations  is  to  make  the  mean 
field  approximation  The  mean  field  approximation  claims  since  we  are  going  to  working 
with  such  large  numbers  of  particles  in  the  ground  state,  we  can  approximate  ar_n  and  at 
as  a  number,  s/Nq.  Normally,  a^=0  would  return  a  value  of  \/N0  —  1.  [17]  The  difference 
between  a  single  particle  in  a  state  consisting  of  100,000  is  insignificant.  When  our  quantum 
numbers  p,  q,  and  k  equal  zero,  all  the  operators  go  to  This  leads  to  a  factor  of 

appearing  in  front  of  the  interaction  term.  NqV(0)/2Ld  is  the  coefficient  in  front  of  the 
zeroth  order  interaction  portion  of  the  Hamiltonian.  This  component  of  the  Hamiltonian 
models  the  interactions  between  the  condensate  atoms  and  other  condensate  atoms.  Next 
we  should  set  three  operators  to  a0.  When  these  three  operators  all  go  to  \/Ao,  we  will 
only  have  one  operator  left.  However  the  expectation  value  of  a  single  operator  is  always 
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zero.  This  stems  from  the  fact  that  an  expectation  value  needs  two  equal  energy  states 
with  the  same  particle  number  or  the  expectation  value  goes  to  zero.  When  a  raising  or 
lowering  operator  acts  upon  our  initial  state,  particle  number  will  change.  Since  we  will 
have  same  momentum  states  with  with  different  particle  number  ,  the  expectation  value  is 
zero.  As  a  result,  H3  will  also  go  to  zero.  Since  we  need  to  return  Fock  states  with  the  same 
particle  number  in  order  to  calculate  with  them,  no  system  of  2  raising  and  1  lowering  or 
2  lowering  and  1  raising  operator  will  return  our  initial  particle.  Hi  =  H3  =  0.  H2,  which 
models  the  interactions  between  the  condensate  atoms  and  the  excited  gas,  will  be  non-zero. 
We  send  two  operators  from  a0  to  \/iVo,  and  keep  two  to  act  upon  the  state.  There  are 
six  combinations  two  operator  systems  when  two  operators  are  approximated  to  ao-  The 
first  term  is  when  the  two  creation  operators  are  set  to  zero,  q  =  p  and  k  =  —p.  This 
leaves  us  with  N0/2LDV(p)a^pap.  The  N0  is  due  to  the  two  ao  we  approximated,  and  the 
remaining  operators  are  summed  over  one  quantum  number.  A  full  derivation  can  be  found 
in  Appendix  A.  After  completing  the  decomposition  for  all  six  possible  combinations,  we 
find 


®  +  +  24«fc)  +  (14) 


k 


E.  Hartree-Fock  Decomposition 


H4  which  describes  the  thermal  cloud  interacting  with  other  thermal  cloud  atoms.  This 
portion  of  the  Hamiltonian  describes  our  momentum  numbers,  k ,  q ,  and  p,  when  none  can 
equal  zero.  Only  atoms  in  excited  states  contribute  to  Hi.  When  we  add  Ho ,  H2,  and  Hi,  we 
have  the  Bogoliubov  Hamiltonian.  This  Hamiltonian  can  describe  an  interacting  condensate 
and  its  surrounding  gas.  The  interaction  component  of  the  Hamiltonian  is  now  mostly 
in  quadratic  operators,  but  it  is  not  yet  diagonal  so  we  cannot  compute  the  eigenvalues 
for  several  reasons.  In  order  to  be  diagonal,  the  Hamiltonian  must  have  no  off-diagonal 
elements.  A  series  of  operators  with  equal  numbers  of  creation  and  annihilation  operators  is 
diagonal  while  having  an  inequality  between  the  number  of  operators  constitutes  being  off 


diagonal. 
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remains  undiagonalized  and  quartic  in  terms  of  its  operators.  We  cannot  use  a  mean 
field  approximation  on  H4,  because  none  of  the  operators  can  go  to  do  because  we’ve  already 
eliminated  these  terms.  Therefore,  we  need  to  use  a  different  method  of  obtaining  an  equiva¬ 
lent  quadratic  form.  Hartree-Fock  decomposition  will  achieve  exactly  this.  I11  Hartree-Fock, 
you  take  every  possible  ordered  pair  of  operators  and  create  expectation  values  of  them, 
leaving  the  other  two  operators  unchanged.  For  our  four  operator  term,  this  will  leave  six 
terms.  Once  we  have  our  six  terms,  we  set  equivalences  that  will  reduce  the  number  of 
quantum  numbers  to  two.  The  ata^  operators  and  our  densities,  normal  and  anomalous, 
stay  in  terms  of  the  quantum  number  k.  p  and  k  are  allowed  to  appear  without  reservation 
in  the  other  portions  of  the  term.  We  define  our  normal  density  to  be  the  expectation  value 
of  n(k )  =  ata£,  because  it  results  in  a  number  operator  and  when  combined  with  1/LD,  will 
provide  the  actual  density  for  our  condensate.  [17]  We  also  define  an  anomalous  density  to 
be  m{k)  =  %%.  After  applying  these  definitions,  our  i74  turns  into 

H4  ~  nV (0)  Y  E  V “  *)  +  rh{q)a +  h.c.)  (15) 

kq^O 

We  now  redehne  our  quadratic  (but  off  diagonal  Hamiltonian)  to  be  H'2  =  H2  +  H 4.  This 
equation  can  be  made  much  simpler  by  substituting  the  self  energies  En(/c),  S20 (k),  and 
Eg2(k).  [18]  We  set  Sn(fe)  equal  to  all  components  of  the  already  diagonalized  portion.  We 
set  S2o (k)  and  its  Hermitian  Sq 2(k)  equal  to  all  the  components  in  front  of  the  off  diagonal 
operators.  This  simplifies  our  equation  and  produces  the  result  below 


k^O 


—  p  +  £n(/c) 


4% 


1 

+  2 


E 


E20(k)aia_-k  +  T.’„2(katat 


(16) 


Where 


£11  (k) 


=  nV  (0)  +  n0V  (k)  +  JqYj  (?  -  k) 

0 

(17) 

£20 (k)  =  n0V (k)  +  Y  m(q)V (?  -  k) 

(18) 

If  we  did  not  want  to  include  H±  in  our  calculations,  we  could  obtain  the  results  from 
only  H2  by  setting  our  densities,  n(q)  and  fh(q),  to  zero  and  moving  11  to  no-  This  result 


20 


would  account  for  interactions  between  the  condensate  and  itself  and  the  condensate  and 
the  thermal  gas  surrounding  the  condensate  (but  no  interactions  between  the  gas  and  other 
thermal  atoms). 


F.  Bogoliubov  Transformation 

Next  we  diagonalize  H2  in  order  to  extract  its  eigenvalues.  Diagonalized  matrices  are 
necessary  for  calculations,  because  we  need  to  find  the  eigenstates  of  the  Hamiltonians  to 
find  the  energies  of  the  system  to  plug  into  the  Bose-Einstein  distribution  (Eq  3).  [19]  An 
off  diagonal  Hamiltonian  has  off  center  elements  such  as  quadratic  operators  which  feature 
two  creation  operators.  Diagonal  elements  need  equal  numbers  of  creation  and  annihilation 
operators.  We  use  the  Bogoliubov  transformation  in  order  to  achieve  this.  The  Bogoliubov 
transformation  uses  a  change  of  variable  to  change  any  non-diagonal  quadratic  equation  into 
a  diagonal  one.  [15]  Unfortunately,  the  transformation  does  not  conserve  particle  number  so 
we  introduce  the  grand-canonical  ensemble  and  n  to  conserve  particle  number.  Bogoliubov 
originally  invented  this  method  in  order  to  solve  mathematical  difficulties  with  supercon¬ 
ductivity.  In  order  to  do  this,  Bogoliubov  redefined  the  creation  and  annihilation  operators 
as  follows 


k  =  ukk  +  v*jP_%  (19) 

&l  =  ulPk  +  v^-k  (2°) 

Instead  of  our  previous  operators  that  create  or  destroy  single  particles,  our  new  operators 
feg  and  bt  create  and  annihilate  quasi-particles.  Quasi-particles  are  a  mathematical  tool  used 
to  describe  complicated  systems.  Instead  of  modeling  a  difficult  system  of  many  moving 
particles  strongly  interacting  with  each  other,  we  can  model  quasi-particles  which  act  more 
like  non-interacting  particles.  Keep  in  mind  quasi-particles  are  not  ’’real”  particles  and 
instead  a  mathematical  representation  used  to  simplify  the  system. 

To  ensure  our  operators  commute,  we  calculate  the  commutation  relation  we  previously 
defined  for  a %  and  at.  As  opposed  to  a  simple  Dirac  delta,  we  find  the  result  u^ut  —  v-v^  =  1. 
First,  we  plug  in  our  definitions  for  ag  and  at  into  the  Bogoliubov  Hamiltonian,  H2 .  Next,  we 
bundle  similar  terms;  We  place  all  bj-b £,  bVb^,  and  btbi  into  separate  groups.  We  find  that  the 
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bundle  of  b and  btbt  are  conjugates  of  each  other.  ( v 2  +  u2  +  2 uv)V(q)  +  2V (0 )vu  +  —  n 

lies  in  front  of  b^b^.  We  have  the  freedom  to  choose  u  and  v  such  that  these  off  diagonal 
elements  vanish.  We  linearly  combine  our  b^b^s  coefficients  with  the  boson  commutation 
relation  via  Mathematica  in  order  to  solve  for  u  and  v.  [4] 


uk  = 


\ 


1  6t 


Vk  = 


A 


1/67 


/i  +  Sn(fc)  +  A 

w k  ) 

h  +  ^ll(fc)  _  A 
w k  j 


where 


^k  = 


We  redehne  our  normal  and  anamolous  energies  using  our  u  and  v. 


(21) 


(22) 


(23) 


n(k)  =  ( u\  +  vl)  f%  +  v\  (24) 

m(k)  =  u%v £  (2  f%  +  1)  (25) 

After  plugging  our  values  for  u  and  v  into  the  our  Hamiltonian,  simplihcation  brings 
us  the  full  diagonalized  Hamiltonian.  Despite  the  complexity  of  our  system,  H'2  ultimately 
arrives  at  a  simple  and  elegant  final  form. 


«2  =  E'0  +  J2^%  (26) 

k+ 0 

G.  The  Hugenholtz-Pines  Theorem 

We  have  yet  to  find  an  expression  for  the  chemical  potential,  ji.  We  can  find  the  chem¬ 
ical  potential  by  introducing  the  mean  held  expansion  into  the  equation  of  motion  for  an 
operator,  a % 
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ihdtdj: 


d^H 

(efc  -  »)  H  +  Jo  ^ 


(27) 


H  —  nV(Q)  +  y~]  y  (fc)  ^n(fc)  +  m(fc))  (28) 

k 

In  1959,  by  a  rigorous  and  general  argument,  Hugenholtz  and  Pines  determined  that  /i, 
the  chemical  potential,  must  equal  the  difference  between  the  two  self  energies  describing 
the  bosonic  system;  [20] 


/i  —  En(0)  —  E20(0)  (29) 

is  the  quasi-particle  spectrum  that  describes  the  dispersion  of  non-condensate  atoms. 
At  vanishingly  small  k,  this  spectrum  must  produce  a  gapless  Goldstone  mode,  meaning 
ojq  =  0.  This  is  due  to  the  continuous  symmetry  of  the  Hamiltonian  that  is  broken  by  the 
presence  of  a  condensate.  One  can  easily  check  that  for  to  be  gapless  at  k  —  0,  the  above 
equation  must  be  satisfied. 


H.  The  Popov  Approximation 

We  can  check  that  when  we  set  k  =  0  in  £n(0)  and  £2q(0),  we  get 


=  nV( 0)  +  n0V( 0)  +  ^  n{q)V(q) 

g^O 

(30) 

S2o(0)  =  n0V (0)  +  ^(9)^(9) 

&  0 

(31) 

We  substitute  the  values  gained  from  Eq.  30  into  the  Hugenholtz-Pines  theorem  (Eq. 
29),  we  arrive  at  the  result 

£n(0)  -  S20(0)  =nV(  0)  +  ^D^V(k)  (n{h)  -  m(fc)) 

k 

7^' 


(32) 
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We  remember  our  value  for  y  in  Eq.  28.  Unfortunately,  our  equations  don’t  immediately 
satisfy  the  Hugenlioltz-Pines  theorem.  This  is  a  problem,  because  it  indicates  our  method 
is  somehow  flawed.  However,  we  can  implement  the  common  Popov  approximation  in  order 
to  comply  with  Hugenholtz-Pines.  The  Popov  approximation  amounts  to  neglecting  the 
anomalous  density,  setting  rh(k )  — >  0.  [4,  17] 

Physically,  rh  corresponds  to  two  non-BEC  atoms  simultaneously  scattering  out  of  the 
BEC  or  scattering  into  it.  This  interaction  happens,  but  is  not  common  enough  to  signih- 
cantly  alter  our  model  of  the  BEC. 

I.  Two  Dimensions:  the  Quasi-Condensate 

As  we  saw  in  section  If  B  2,  when  we  consider  a  gas  of  polarized  dipolar  bosons,  their 
interactions  are  anisotropic.  When  the  dipoles  are  polarized  along  the  z-axis,  the  interactions 
are  repulsive  in  the  x-y  plane,  but  are  attractive  in  the  ^-direction,  so  two  dipoles  will  attract 
if  aligned  in  a  “head-to-tail”  configuration.  For  this  reason  it  is  advantageous  to  confine  the 
dipolar  atoms  in  the  x-y  plane  with  some  sort  of  strong  trapping  potential,  like  a  harmonic 
trap  in  the  ^-direction,  to  discourage  the  gas  from  collapsing  onto  itself.  Experimentally, 
this  trapping  potential  can  be  created  using,  for  example,  strong  laser  beams  that  attract 
the  atoms  into  a  plane.  This  results  in  an  effectively  two  dimensional  system  where  all  of  the 
atoms  occupy  a  wave  function  with  a  common  z-dependence  given  by  y(z),  which  we  take 
to  be  a  Gaussian  (the  single  particle  ground  state  of  a  ID  harmonic  potential).  We  discuss 
this  further,  and  derive  the  reduced  contact  and  dipolar  interactions  that  are  relevant  for  a 
quasi-2D  geometry  in  Appendix  E. 

When  Bose  gas  is  forced  into  2D,  it  turns  out  that  a  true  BEC,  where  every  atom  occupies 
the  exact  same  quantum  mechanical  ground  state,  can  only  exist  at  exactly  zero  temperature. 
Thus,  no  BEC  can  exist  in  2D  at  any  finite  temperature  [21],  Instead,  however,  a  different 
low  temperature  phase  can  emerge,  called  a  “quasi-condensate,”  which  looks  like  a  BEC  on 
small  length  scales,  but  has  large  thermal  fluctuations  in  the  phase  of  the  condensate  wave 
function  on  longer  length  scales.  This  fact  was  proven  rigorously  by  Mermin  and  Wagner,  and 
is  thus  called  the  Mermin- Wagner  theorem.  [22]  It  states  that  in  2D,  continuous  symmetries 
cannot  be  spontaneously  broken  at  finite  temperature  in  systems  with  sufficiently  short- 
range  interactions  (which  includes  both  contact  and  dipolar  oc  1/r3  interactions).  In  this 
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case,  the  continuous  symmetry  is  the  phase  symmetry;  in  our  Hamiltonian,  the  operators 
can  be  modihed  by  any  phase  e1^  and  the  Hamiltonian  remains  unchanged.  Essentially,  any 
long-range  interaction  can  induce  fluctuations  with  minimal  energy  costs,  hence  altering  the 
ground  state  of  the  system. 

Heuristically,  we  find  that  the  Hartree  Fock  Bogoliubov  (HFB)  equations  can  not  be 
solved  at  finite  temperatures  in  2D,  due  to  the  divergence  of  the  integrals  over  h{k).  This  is 
a  backwards  proof  of  the  Mermin- Wagner  theorem.  So,  we  need  to  introduce  a  new  mean- 
field  for  the  ground  state  phase.  Instead  of  a  pure  BEC,  we  will  have  a  quasi-condensate. 
The  BEC  was  described  by  the  k  =  0  (spatially  uniform)  mode  ^/uq.  To  determine  an 
analogous  mode  for  the  quasi-condensate,  we  consider  the  number-phase  representation  of 
the  Bose  annihilation  and  creation  operators  that  we  started  with.  We  can  write 


T(r)  =  v^e^.  (33) 

where  n(f )  is  the  number  operator  and  9(r)  is  the  phase  operator.  Anticipating  there  being 
fluctuations  in  the  phase  in  a  quasi-condensate,  we  will  leave  the  phase  operator  as-is,  but 
introduce  a  mean-field  (just  a  number)  for  the  mean  value  of  the  density, 

T(r)  =  y/no  +  5n(r)e1^.  (34) 

Expanding  this  to  linear  order  in  the  operators,  which  we  have  to  do  to  match  up  with 
our  theory  above,  we  find 


'k(r)  =  V^o  +  iyfa6(r)  + 

l  yjTLQ 

itr)  =  ,/SS  -  (35) 

We  now  wish  to  make  a  connection  with  the  Bogoliubov  theory  outlined  above  in  sec¬ 
tion  11 F;  we  move  to  momentum  space,  and  proceed  to  represent  these  operators  in  momen¬ 
tum  space.  Enforcing  the  commutation  relation  for  number  and  phase  operators, 


5n(r'),  9(r) 


U(r'  -  r), 


(36) 


25 


we  find  the  momentum  space  forms  to  be 


%  =  +  vk )  (h  + 


9%  2\/n0 


(' uk  ~  vk )  b 


(37) 


The  normal  density  of  non-condensate  atoms  can  be  written  as 


n(k)  =  n0(%%)  +  -^-{5h%8h%)  +  \ 

=  Mhh)  + 


h 


(38) 


Let’s  look  at  the  behavior  of  this  function  at  small  values  of  k.  We  find  that  it  goes  like 
h{k)  ~  T/k 2  at  small  k.  In  3D,  this  is  convergent  at  small  k  when  integrating  over  n(k). 
In  2D,  however,  this  is  divergent!  We  notice  that  it’s  the  term  in  particular  that 

behaves  this  way.  This  makes  sense,  because  this  term  involves  the  phase  operator,  and 
we  expect  that  phase  fluctuations  to  be  large  and  kill  the  BEC  in  2D.  Thus,  we  define  the 
quasi-condensate  as 


n0^n0  +  n0J2(hh)-  (39) 

k 

Thus,  the  quasi- condensate  absorbs  the  phase  fluctuations  that  kill  the  BEC.  The  quasi¬ 
condensate  still  has  highly  suppressed  density  fluctuations,  however. 

To  fix  the  HFB  equations,  we  can  simply  take  every  term  involving  no(%%)  and  set  it 
to  zero.  This  gives  convergent  integrals,  and  a  well  behaved  theory  for  a  quasi-condensate 
that  is  perfectly  analogous  to  the  3D  theory  for  the  true  BEC. 

J.  Superfluidity 

One  of  the  most  interesting  phenomenon  associated  with  BECs  is  the  emergence  of  su¬ 
perfluidity.  Superfluids  are  fluids  which  experience  non-viscous  flow,  are  interesting  from 
a  practical  standpoint  because  easy  access  to  superfluids  would  eliminate  major  losses  of 
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energy  and  work  for  industrial  applications.  Heuristically,  it  has  been  found  that  the  super¬ 
fluidity  coincides  with  the  existence  of  certain  types  of  BECs.  However,  not  every  BECs  is 
a  superfluid  and  not  all  superfluids  stem  from  the  existence  of  a  BEC. 

The  properties  of  superfluid  BECs  were  being  studied  even  before  the  connection  between 
BECs  and  superfluids  was  known.  Lev  Landau,  1962  Nobel  Prize  winner,  analyzed  Helium 
4,  a  boson,  for  its  superfluid  properties.  He  showed  that  Helium  4  would  not  produce  sound 
wave  excitations  while  flowing  past  a  wall,  granted  the  flow  of  the  Helium  was  slower  than 
the  sound  velocity.  Landau  suggested  the  sound  velocity  was  also  the  critical  velocity  for 
superfluidity.  Above  this  velocity,  Helium  4  is  not  a  superfluid.  Unknown  at  the  time,  this 
phenomenon  was  due  to  the  existence  of  a  BEC  within  Helium  4.  [23] 

In  order  for  a  BEC  to  support  superfluidity,  there  needs  to  be  interactions  between  its 
constituents.  [22]  A  non-interacting  BEC  is  not  a  superfluid.  Different  types  of  interactions, 
whether  contact  or  dipolar,  could  therefore  produce  superfluidity,  and  the  superfluid  fraction 
(the  density  of  particles  acting  as  a  superfluid  over  the  total  density  of  particles)  will  vary 
depending  on  it.  Not  only  do  interactions  factor  into  the  superfluid  (SF)  fraction,  but 
the  degrees  of  freedom  influence  it  as  well.  We  calculate  the  superfluid  fraction  using  the 
following  equation 


(40) 


III.  RESULTS 


A.  Three  Dimensions  with  Contact  Interactions 


As  explained  previously,  a  BEC  with  no  interactions  will  correspondingly  have  no  super- 
fluid  properties.  Therefore,  the  non-interacting  case  in  three  dimensions  shown  in  Fig.  1 
from  the  Non-Interacting  Bose  Gas  Section  shows  only  the  condensate  fraction.  It  has  no 
superfluid  properties  so  the  superfluid  fraction  remains  at  zero. 

Contact  Interactions  bring  significantly  more  interesting  results  than  the  non-interacting 
case.  From  seeing  both  the  condensate  fraction  and  superfluid  fraction  on  a  single  plot, 
we  see  that  contact  interactions,  no  matter  how  minimal,  will  give  the  SF  fraction  a  slight 
increase  over  the  BEC  fraction.  As  we  increase  the  interaction  strength,  we  similarly  find 
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FIG.  2:  Condensate  fraction  and  superfluid  fraction  of  a  3  Dimensional  BEC  with  contact  inter¬ 
actions  of  increasing  strength.  As  the  contact  interaction  strength  increases,  so  does  the  BEC 
fraction,  superfluid  fraction,  and  critical  temperature.  The  superfluid  and  BEC  have  the  same 
critical  temperature.  As  temperature  nears  zero,  BECs  with  very  strong  contact  interactions  ex¬ 
perience  the  phenomenon  known  as  quantum  depletion  in  which  some  atoms  occupy  the  thermal 
cloud  even  at  zero  0  Kelvin. 

a  greater  increase  in  the  relative  increase  of  the  SF  fraction  over  the  BEC  fraction.  We 
also  increasing  the  contact  interaction  strength  increases  the  critical  temperature.  In  3D 
with  contact  interactions,  the  superfluid  and  condensate  crtical  temperatures  are  the  same. 
Near  absolute  zero,  we  see  the  effect  of  quantum  depletion  diminishing  the  condensate. 
Quantum  depletion  occurs  because  of  the  large  strength  of  the  interactions,  particles  with 
non-zero  momenta  are  pushed  out  of  the  condensate,  even  at  zero  temperature.  The  system  is 
stabilized  with  excited  atoms  at  zero  temperature.  On  the  other  end  of  the  temperature  scale, 
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FIG.  3:  Condensate  fraction  and  superfluid  fraction  in  a  nearly  2D  quasi-condensate  at  finite  tem¬ 
peratures  with  varying  strengths  of  contact  interactions.  As  contact  interaction  strength  increases, 
so  do  both  the  superfluid  and  condensate  critical  temperatures:  the  superfluid  critical  tempera¬ 
ture  remains  the  lower  of  the  two.  The  occupancy  of  the  condensate  and  superfluid  also  increase 
dramatically  with  increasing  contact  interaction  strength. 

we  find  that  the  critical  temperature  is  extended  significantly  beyond  the  non-interacting 
critical  temperature.  In  fact,  every  increase  in  contact  interaction  strength  will  induce  a 
slightly  larger  critical  temperature. 


B.  Two  Dimensions  with  Contact  Interactions 


In  2D,  a  pure  BEC  only  exists  at  0  Kelvin.  However,  the  quasi- condensate  functions 
similar  to  a  pure  condensate.  The  quasi-condensate  has  supressed  density  fluctuations  and 
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FIG.  4:  Critical  temperature  of  the  superfluid  fraction  and  the  quasi-condensate  fraction  at  different 
contact  interaction  strengths  while  in  two  dimensions.  The  critical  temperature  for  the  superfluid 
and  the  quasi-condensate  both  initially  sharply  increase  and  then  continue  to  increase  more  slowly 
with  increasing  contact  interaction  strength.  However,  the  superfluid  critical  temperature  spectrum 
remains  moderately  lower  than  the  quasi-condensate  critical  temperature  spectrum. 

rather  large  phase  fluctuations.  Mathematically,  these  phase  fluctuations  manifest  in  the 
term  f  +  1).  This  behaves  like  p  when  taken  at  small  k.  While  this  is 

easily  handled  in  three  dimensions,  this  term  diverges  in  2D.  However,  we  absorb  it  into 
no,  the  condensate  density.  This  converts  our  condensate  from  strictly  a  BEC  to  a  quasi¬ 
condensate  with  phase  fluctuations.  When  taken  from  a  small  distance  away,  a  BEC  and 
quasi-condensate  act  the  same.  However,  at  large  length  scales,  a  quasi-condensate  will 
change  phase  whereas  a  BEC  has  a  single  phase.  For  the  purposes  of  superfluidity  and 
condensate  fraction,  two  dimensions  have  significant  effects. 


30 


Similar  to  3D,  increasing  contact  interaction  strength  will  likewise  increase  the  condensate 
fraction.  The  only  temperature  where  lower  interaction  strength  leads  to  a  higher  condensate 
fraction  is  near  absolute  zero  when  strong  interactions  cause  quantum  depletion.  However, 
superfluid  fraction  is  only  increased  with  greater  interaction  strengths.  The  increase  in 
contact  interaction  strength  plays  a  significantly  greater  role  in  the  increase  of  the  condensate 
and  SF  fraction  than  in  3D.  Unlike  in  3D,  the  superfluid  has  a  different  critical  temperature 
than  the  condensate.  The  condensate  undergoes  a  superfluid  transition  initially  and  enters 
a  phase  where  the  condensate  still  exists,  but  does  not  possess  superfluid  properties.  Then 
the  condensate  undergoes  a  transition  at  its  critical  temperature  into  a  thermal  gas.  The 
two  transitions  can  be  seen  explicitly  in  Fig.  4. 

A  derivation  of  the  quasi-2D  contact  interaction  potentials  is  given  in  Appendix  D. 

C.  Two  Dimensions  with  Dipolar  Interactiions 

The  dipolar  interaction  affects  the  condensate  significantly  more  dramatically  than  con¬ 
tact  interactions.  Although  weakly  interacting  dipoles  provide  similar  condensate  fractions 
and  excitation  spectrum  as  those  of  contact  interactions,  strong  dipolar  interactions  cause 
large  changes  in  the  condensate’s  nature.  Similar  to  2D  contact  interactions,  an  increase 
in  the  dipolar  interaction  strength  produces  an  increase  in  the  condensate  critical  tempera¬ 
ture.  However,  strong  dipolar  interactions,  where  rotons  are  present,  such  as  the  green  line 
in  Fig.  6  show  an  enormous  quantum  depletion.  For  a  moderately  strong  interaction,  l  of 
the  total  particles  are  in  the  thermal  cloud  at  zero  temperature!  Weak  Dipolar  interactions, 
where  rotons  do  not  manifest,  produce  higher  condensate  fractions  over  the  majority  of  the 
temperature  scale. 

Superfluids  also  undergo  an  interesting  transition.  While  weakly  interacting  dipoles  pro¬ 
duce  an  increase  in  superfluid  fraction  and  superfluid  critical  temperature,  the  development 
of  the  roton  excitation  spells  disater  for  superfluidity.  Not  only  does  superfluidity  decrease 
dramatically,  but  the  roton  excitation  causes  the  lower  superfluid  fraction  than  more  weakly 
interacting  dipoles.  In  fact,  the  spectrum  of  superfluid  fraction  for  strong  dipoles  concaves. 
Strong  dipolar  interactions  and  roton  excitations  lead  to  a  dramatic  decrease  in  the  conden¬ 
sate’s  ability  to  support  superfluidity. 

Rotons  have  a  profound  effect  on  quantum  depletion.  As  the  condensate  experiences  a 
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FIG.  5:  Effect  of  dipolar  interaction  strength  on  condensate  occupancy  at  0  Kelvin.  The  condensate 
experiences  quantum  depletion  at  a  consistent  rate  until  dipolar  interactions  become  strong,  at 
which  point  the  roton  mode  appears  and  cause  a  rapid  depletion  in  the  condensate. 

gradual  increase  in  dipolar  interaction  strength,  a  slight  linear  drop  in  condensate  fraction 
occurs.  Near  the  dipolar  interaction  strength  of  gSHD /V2nlz  =  1.5,  the  roton  mode  begins 
to  form  and  creates  a  sharp  decrease  in  the  condensate  occupation  at  0  Kelvin.  This  effect 
is  clearly  illustrated  in  Fig.  5.  This  steep  drop  continues  until  near  gd^D/V^h  =  3 
at  which  point  the  condensate  fraction  at  0  Kelvin  begins  to  return  to  smaller  decreases 
with  increasing  interaction  strength.  We  can  see  the  effect  of  dipolar  interaction  strength  in 
regards  to  excitation  energy  out  in  Fig.  7.  This  probably  occurs  because  a  roton  mode  forms 
near  gd^iol'J 2ttIz  =  1-5  which  causes  lots  of  low  energy  excitations.  This  process  compounds 
with  increasing  dipolar  interaction  strength  and  continues  to  shrink  the  condensate  fraction 
at  zero  temperature. 
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The  spectrum  for  strong  dipolar  excitations  is  unique  because  with  increasing  momentum, 
we  find  a  local  maximum  (a  small  hill),  and  then  a  wide  minimum  where  the  roton  excitation 
exists.  If  we  have  only  weakly  interacting  dipoles,  then  the  excitation  spectrum  looks  similar 
to  those  of  contact  interactions.  Fig.  8  shows  that  only  at  low  temperatures  do  roton 
modes  exist.  As  temperature  increases,  many  atoms  become  excited  out  of  the  condensate, 
stabilizing  it.  If  we  calculated  without  a  self  consistent  solution,  this  would  not  occur  and 


FIG.  6:  Condensate  fraction  and  superfluid  fraction  of  a  2  Dimensional  BEC  with  dipolar  inter¬ 
actions  of  varying  strength.  Increasing  dipolar  interaction  strength  vastly  increases  the  critical 
temperature  of  the  condensate.  However,  it  also  drastically  increases  the  effect  of  quantum  de¬ 
pletion,  causing  condensates  with  weak  dipolar  interactions  to  have  a  greater  condensate  fraction 
at  most  temperatures.  Increasing  dipolar  interaction  strength  initially  has  a  positive  influence  on 
the  superfluid  critical  temperature  and  superfluid  fraction,  but  too  large  of  dipolar  interactions 
decreases  the  critical  temperature  and  superfluid  fraction. 
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FIG.  7:  Comparison  of  self  consistent  versus  non-self  consistent  method  of  calculating  in  regards 
to  excitation  energy  of  increasing  dipolar  interaction  strength  at  0  K.  When  calculating  not  self 
consistently,  the  roton  excitation  energy  drops  to  zero  soon  after  roton  excitations  appear.  There¬ 
fore,  the  condensate  collapses  and  all  atoms  spill  into  the  thermal  cloud.  However,  a  self  consistent 
calculation  shows  that  regardless  of  the  interaction  strength,  the  roton  excitation  energy  will  never 
go  to  zero.  As  a  result,  the  system  is  mechanically  stable. 

instead  the  condensate  would  collapse  due  to  the  strong  dipolar  interactions,  regardless  of 
temperature.  However,  we  find  that  temperature  acts  as  a  stabilizing  mechanism  for  the 
condensate,  removing  the  roton  mode. 

A  derivation  of  the  quasi-2D  dipolar  interaction  potentials  is  given  in  Appendix  E 
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FIG.  8:  Excitation  spectrum  of  dipolar  interactions  at  varying  temperatures  and  occupation  of  the 
thermal  cloud  at  varying  temperatures  and  increasing  momentum.  As  temperatures  increase,  the 
roton  excitation  requires  more  energy.  With  enough  temperature,  the  roton  excitation  disappears 
entirely.  The  thermal  cloud  occupation  reveals  this  feature  by  revealing  a  larger  density  of  atoms 
at  the  roton  momentum  at  zero  temperature.  As  the  temperature  rises,  the  overall  density  of 
excitations  increases  as  well  as  the  majority  of  the  excitations  moves  away  from  the  former  roton 
minimum  towards  lower  momentum  excitations. 

IV.  CONCLUSION 

Particles  at  ultracold  temperatures  have  two  important  features,  their  wavelike  nature 
and  whether  the  particles  are  bosons  or  fermions.  Bosons,  the  particles  we  model,  are  integer 
spin  particles  that  can  occupy  the  same  quantum  space.  Fermions,  particles  of  half  integer 
spin,  cannot.  We  model  boson-like  atoms  which  consist  of  an  even  number  of  fermions,  such 
as  Sodium  23.  As  we  cool  a  gas  of  bosons  into  the  ultracold  regime,  we  find  the  De  Broglie 
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wavelength  of  the  atoms  begins  to  increase.  Once  the  size  of  the  wavelength  becomes  on  the 
order  of  the  inter  particle  spacing,  a  Bose  Einstein  condensate  forms.  This  unique  statistical 
phenomenon  creates  a  new  state  of  matter  consisting  of  a  macroscopic  quantum  mechanical 
state.  Once  the  gas  is  cooled  to  absolute  zero,  every  particle  enters  the  BEC  (only  in  the 
non-interacting  case!). 

Specifically,  we  modeled  dipolar  interactions  within  the  BEC.  However,  dipolar  interac¬ 
tion’s  attractive  potential  causes  the  condensate  to  collapse  when  in  three  dimensions  so  we 
confine  our  BEC  to  a  nearly  two  dimensional  geometry.  As  a  result,  the  condensate  gains 
moderate  phase  fluctuations,  creating  a  quasi-condensate.  We  treat  the  quasi-condensate 
with  nearly  identical  mathematical  formalism,  but  still  accounting  for  the  phase  fluctua¬ 
tions.  The  excitation  spectrum  for  dipolar  interactions  in  nearly  two  dimensions  creates  a 
peculiar  excitation  called  a  roton.  The  roton  allows  for  high  momentum  excitations  for  very 
low  energy.  We  investigate  the  extent  to  which  rotons  influence  our  condensate  system. 

Due  to  the  Schrodinger’s  equation  restrictions  on  symmetric  wave  functions  and  degrees 
of  freedom,  we  find  it  prohibitively  difficult  to  solve  this  problem.  Therefore,  we  use  an 
alternative  method  known  as  2nd  Quantization.  We  enumerate  number  states  and  apply 
creation  and  annihilation  operators  to  them  in  order  to  count  how  many  particles  we  have 
in  each  Fock  state.  We  further  make  the  mean  held  approximation,  Hartree-Fock  decompo¬ 
sition,  Bogoliubov  transformation,  and  finally  the  Popov  approximation  in  order  to  create 
a  diagonalized  quadratic  Hamiltonian. 

We  fold  that  in  nearly  two  dimensions,  strong  dipolar  interactions  cause  a  significant 
decrease  in  the  superfluid  fraction.  In  addition,  two  dimensions  cause  two  phase  transitions, 
first  from  a  superfluid  condensate  to  a  non-superhuid  condensate  and  second  from  a  non- 
superhuid  condensate  to  a  thermal  gas.  When  we  calculate  our  solution  self  consistently, 
ensuring  particle  number  stays  constant,  we  find  that  the  condensate  will  not  collapse  at 
zero  temperature  regardless  of  the  dipolar  interaction  strength.  While  non-self  consistent 
calculations  cause  the  roton  excitation  energy  to  go  to  zero,  exciting  every  atom  out  of 
the  condensate,  self  consistent  calculations  reveal  the  roton  energy  never  goes  to  zero.  The 
system  remains  mechanically  stable.  In  addition,  raising  the  temperature  of  our  system 
gradually  removes  the  roton  excitation,  stabilizing  the  system.  Our  results  will  be  revealing 
for  modern  experiments  with  ultra  cold  atoms  focusing  on  dipolar  interactions  and  the 
condensate. 
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Appendix  A:  Derivation  of  Hartree-Fock  Equations 

We  begin  with  the  total  Hamiltonian  describing  our  system.  We  have  our  non-interacting 
component  as  the  first  term,  which  describes  the  collective  kinetic  energy  of  the  system.  Our 
second  term  consists  of  the  potential  term  which  will  vary  depending  on  which  potential  we 
introduce. 


H  —  ^\r)Hl^>{f)df  +  -  /  dd(r)dd(r')U(r  —  r')d'(r)d'(r') 


(Al) 


Tis  a  creation  operator.  It  will  create  a  particle  at  a  certain  position  of  r.  Using  an 
orthonormal  basis  of  plane  waves,  we  set  the  creation  operator  T  to  another  operator  a £ 
which  creates  a  particle  of  a  certain  momentum.  An  orthonormal  basis  is  a  basis  of  unit 
vectors  which  are  all  orthogonal  to  each  other.  An  example  would  be  x ,  y ,  and  £.  They  are 
three  unit  vectors  which  are  orthogonal  or  perpendicular  to  each  other. 


T(rO  =  V-L 


(A2) 


where  ag  acts  as  an  annihilation  operator,  at  acts  as  a  creation  operator. 
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After  substituting  T  and  dA  into  the  kinetic  energy  portion  of  the  Hamiltonian,  we  find 


kq 


-ik-r 


-Wr  df 


(A4) 


The  Laplace  operator  V2  will  return  a  value  of  —  q2.  After  moving  this  value  outside  the 
integral,  our  equation  simplifies  to 


kq 


2atag- 


e-ik-reiq-rdf 


(AS) 
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Because  of  our  orthonormal  basis, the  k  and  q  vectors  are  perpendicular  and  make  the 
above  integral  zero,  except  when  k  —  q.  When  both  momentums  are  equal,  the  integral 
instead  goes  to  Ld. 


=  k2~a\hk  (A6) 

k 

We  follow  the  same  process  for  the  interaction  potential,  with 
V(r  —  r')  =  EpZd^(p)e~^^_r^  By  following  the  same  process  above,  we 

find 

Hint  =  2^5  ^  V (AT) 

kqp 

At  ultracold  temperatures,  bosonic  gas  forms  a  Bose  Einstein  Condensate,  where  the 
ground  single  particle  state  becomes  macroscopically  occupied  by  N0  atoms.  The  result  of 
ground  state  operators  (k  =  0)  on  our  Fock  states  (|Ao...)  is  approximately  equal  to  y/No 
times  our  original  Fock  state.  Thus,  we  can  introduce  the  mean  field  approximation  found 
below. 


(a0)  ~  (a£)  =  y/iVo  (A8) 

We  introduce  this  approximation  to  our  Hamiltonian  in  A5  and  A6.  We  expand  the 
Hamiltonian  in  the  form  H  =  H0  +  H1+H2  +  ...  where  each  term  has  0,  1,  or  2  non-condensate 
(non-ground  state)  operators.  Noting  that  k  =  0,  and  introducing  the  condensate  density 
of  n0  =  N0/Ld,  we  find 


H0  =  E0=  ^n0V0.  (A9) 

This  term  is  the  condensate  energy  because  it  describes  all  energy  within  the  conden¬ 
sate.  The  next  term  H 1  is  neglected  because  terms  with  1  operator  will  vanish  when  their 
expectation  value  is  taken.  Next  we  move  to  work  on  H2.  There  are  six  possible  ways  to 
set  two  operators  to  k  —  0.  As  an  example,  the  first  two  operators  will  be  set  to  zero  to 
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show  the  process  for  the  mean  field  approximation.  We  set  q  — >  —p  and  k  -+  p.  Our  first 
two  operators  now  equal  a0  =  y/No  and  can  be  brought  out  front.  Our  resulting  term  is 
0%  JW  V(p)  After  making  a  shift  in  quantum  numbers  (p  — >  —k),  our  new  term 

equals 


2^5  Z)  V$)  •  (A10) 

k 

By  setting  each  possible  combination  of  operators  to  zero  and  forcing  our  equations  to 
use  similar  momentum  numbers,  we  find  the  result 


H2 


fa  -  a4)  +  y^[2V +  v<$)  +  H^-k  +  4aw) 


(All) 


The  Hamiltonian  above  does  not  conserve  particle  number  so  we  must  introduce  a  La¬ 
grange  multiplier,  the  chemical  potential,  to  keep  particle  number  constrained.  We  place 
the  chemical  potential,  p  as  a  constraint  on  e. 

The  next  term  in  our  Hamiltonian  expansion  is  H3,  but  just  at  II 1  went  to  zero,  //:*,  also 
will.  Odd  numbers  of  raising  and  lowering  operators,  regardless  of  order  or  capacity,  will 
vanish  when  an  expectation  value  is  calculated.  The  final  component  of  the  Hamiltonian, 
H4:  is  quartic.  At  this  point,  the  rest  of  the  Hamiltonian  is  quadratic  so  we  need  to  force 
H4  to  be  quadratic  as  well.  We  use  Hartree-Fock  decomposition  on  //4  to  accomplish  this. 
We  take  each  combination  of  2  terms  and  take  an  expectation  value  of  them. 
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We  introduce  two  densities  into  this  equation  now.  The  first  is  called  the  normal  density 
hk  and  the  second  is  called  the  anomalous  density  fhk ■ 
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nk  =  (aia%)  (A14) 

mk  =  (a^a_^)  =  (ataf_-)  (A15) 

In  order  to  substitute  the  densities  into  H$,  we  create  equations  to  swap  our  quantum 
numbers.  For  our  first  term  within  the  bracket,  we  set  p  — >  —  k  resulting  in  at  -at  JU^-rhu. 

q — k  k — q 

Following,  we  set  p  — >  k  —  q.  This  action  results  in  atal In  order  to  get  this  term 
into  a  common  form  with  the  other  five  terms,  we  swap  the  quantum  momentum  numbers 
p  -H-  k.  Finally,  our  first  term  is  in  the  desired  form  below. 

Term  1  =  ata*  ^  2LD^ (A16) 

Term  2  will  use  a  normal  density  as  oppose  to  the  anomalous  density  used  in  term  1.  In 
order  to  alter  term  2  into  the  desired  form,  we  set  q  =  0  resulting  in  dtd^(dtd^)V (0) .  After 
substituting  the  normal  density  and  allowing  for  a  quantum  number  swap  p  -B-  k. 

We  follow  this  same  procedure  for  the  remaining  four  terms.  All  process  similar  to  those 
described  above.  Our  new  hf4  follows  below 


1 
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A  =  Efe  -  + 

kyt  0  k^O 

We  take  H2  and  add  it  to  to  form  the  Hamiltonian  for  non-condensate  atoms  inter¬ 
acting  with  non-condensate  atoms. 


H2  -  -  n  +  Tu(k))ald^  +  -  ^  ( T20(k)d^d_ £  +  T*02(k)ald^^j  .  (A19) 


fcyo 


k^  0 


We  find  the  above  equation  by  selectively  defining  En  and  E2q- 
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£n (k)  =  nV(0)  +  -j-^n{q)V{q-  k )  +  n0(0)  +  nQV(k ) 


(A20) 


After  combining  terms,  we  find 


Eli (k)  =  nV (0)  +  n„V(k)  +  k) 


(A21) 


q^O 


Similarly,  we  selectively  choose  £2o  =  £02  to  equal  the  terms  associated  with  aga_g  and 
ata^  After  minimal  simplification, 

k  —k 


£20 (fc)  =  noV (k)  +  ^  m(q)V(q-  k) 

q^O 


(A22) 
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Appendix  B:  Expansion  via  the  Bogoliubov  Transformation 

In  order  to  diagonalize  our  Hamiltonian,  we  use  the  Bogoliubov  Transformation.  Bogoli- 
ubov’s  process  redefines  our  creation  and  annihilation  as  follows 


■nk  +  viP-n 

(Bl) 

?n  * 

?n-*- 

+ 

cr> 

1 

(B2) 

Using  these  relations,  we  also  redefine  our  commutation  relation  in  new  raising  and  low¬ 
ering  operators. 


“it) 


=  8, 


kk1 


(B3) 


ukk  +  u*pbl  +  v%b_% 


=  8 


kk1 


(B4) 


ukkiu*k'b\, 

+ 

ukk’v«b~k' 

+ 

v%]  -,  utjl 
k  -AU  kf  k' 

+ 

t3e 

1 

* 

—  8kk'  (B5) 


By  definition, 


ukhkivk’h-k' 


v*p  tjUpbL 

k  —  AD  Kf  k' 


=  0. 


(B6) 


ukuk' 


Mjpl  +vlvk'  I# 


vkvk' 


-fc’  ~k' 


=  8. 


kk' 


(B7) 


V^Vpdkk'  8kk' 


(B8) 


We  set  k  =  k' ,  forcing  the  Dirac  delta  functions  equal  1.  This  results  in  a  normalization 


constraint  for  u £  and  v%. 


Ve  ~  viv(  =  1 


(B9) 
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Now  we  substitute  our  redefined  creation  and  annihilation  operators  into  our  Hamiltonian. 
For  ease  of  explanation,  we  will  only  show  H2  with  hastened  simplifications. 


#0,2  = 


^  °2Ld^  ^  +  ^2  S  [V$)  («g«!rfc  +  H^-k  +  2S£<%)  +  2V(°)^k 


(BIO) 


#0,2  = 


(iV0^H(0)  +  ^  |V(£)  ((u£6t  +  +  (ujfo  + 


+  2(u*p,  +  v~b_%){u%b%  +  v#l-))  +  2V  (0){u}bi  +  v%b_%){u%b%  +  v}P_t) 


fc  fc  1  uk'J-k>\“1k'Jk  1 


(BH) 


We  multiply  out  all  raising  and  lowering  operators,  combining  like  terms.  We  remove  the 
off  diagonal  terms  and  set  them  equal  to  zero.  We  also  assume  all  our  values  for  k  are  real 
causing  —k  =  k 


2  (ukuk  +  vkvk  +  2ukvk)V(k)  +  (2W)^  (0)  +  u%v-k 


p) 


bkbk  =  0 


(B12) 


The  term  for  b*-:b*.  is  simply  the  complex  conjugate  of  the  above  term.  However,  since 
both  terms  are  real,  the  resulting  term  will  only  be  a  conjugate. 

We  need  to  add  the  single  particle  energy  term,  ek,  into  our  Hamiltonian  to  ensure  we 
have  the  full  spectrum  of  energy  (interaction  potential  and  kinetic  energy).  By  combining 
the  commutation  relation  u2-  —  v2  =  1  and  the  off  diagonal  terms,  we  can  solve  for  the  values 

k  k 

of  u-f, :  and  uy  We  utilize  Mathematica  to  solve  and  partially  simplify.  Through  additional 
algebra  steps,  we  find 


where 


uk 


1  L/i  I  + 

V2  V  w, - 


(B13) 


^  =  4(2  V(k)  +  ec)i 


(B14) 
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=  —sgnV(k)—= 


1  V{k)  +  e^-g 


V2  V  ^ g 

The  remaining  diagonal  terms  are  show  in  their  entirety  below 


(B15) 


Ho, 2  = 


(N0)2V (0)  ,  N0  ^  r 


2LD 


+ 


[v (*)  \ulv^[h  +  vku*kki> l  +  ukv*Ml  +  vluki>lh  +  2uluk% 


+ 


2vkv*kh^l )  +  2V(°)  (uluAk  +  vkvlh^l)  +  ( uluk  +  vkvk*(ek  -  »)) 


k  k  k 


k  k  k 


(B16) 


We  replace  with  bti^  +  1  as  allowed  by  the  commutation  relation.  In  addition,  we 
remove  all  *  because  we  assume  u £  and  v £  are  real.  After  algebraic  simplification  we  find 


H 


0,2 


(JV2^(0)  +  W  £  N>  (“I  +  2uevi  +  “D  &&  +  2V(0)(«|  +  «I)6t6j+ 


^ (fc) (2 upj%  +  vl)  +  2V (0)u|(«|  +  n|) (eg  -  -  //) 


1  ^k^k^k 


(B17) 


We  proceed  to  plug  in  our  values  for  u %  and  v%.  After  simplifying  using  Mathematica, 
our  Hamiltonian  for  Hq^  simplifies  to 


A, = + V,  £  +  ¥  £  h  -  (<*  -  d  -  -n*> 


(Bis) 
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Appendix  C:  Number  Operator  in  terms  of  Creation  and  Annihilation  Operators 

In  our  Hamiltonian,  we  find  pairs  of  quadratic  creation  and  annihilation  operators  in  the 
form  atag.  This  acts  as  a  number  operator  N,  an  operator  which  counts  the  number  of 
particles  in  a  given  state. 

In  order  to  prove  ata^  acts  as  a  number  operator,  we  cat  on  the  following  system 
\Eq,  Ei,  E2l  E3,  E.i, ...  >. 


Nk=2\k0,  h,  k2,  k3,  k4, ...  >=  ki,k2,  k3,  k4, ...  >  (Cl) 

>=  \Zk2aljk0,k1,k2-l,k3,k4:...  >  (C2) 

Vfaatjko,  h,  k2,  k3,  k4, ...  >=  \j k~2\Jk2  -  1  +  l\k0,ku  k2,  k3,  kA, ...  >  (C3) 

Nk=2  =  k2\k0,k1}k2,k3,k4, ...  >  (C4) 


We  can  apply  this  knowledge  to  examining  a  non-interacting  group  of  particles.  We  have 
a  system  of  four  particles 


Xk  =  \k0,  h,  k2,  k3,  k4  >=  |0, 1,  0,  2, 1  > 


(C5) 


When  we  determine  the  energy  for  this  system,  we  apply  the  Hamiltonian  for  a  non¬ 
interacting  system. 


H  =  at  (C6) 

k 

After  counting  each  individual  states  kinetic  energy  and  summing  over  all  the  states,  we 
find  our  collective  energy. 


Hx~k 


EXk 


k\  2  kj  k'j\ 
2  2  2  ) 


10,1,0,2,1  > 


(C7) 
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Appendix  D:  Derivation  of  Quasi-2D  Contact  Interaction  Potential 


We  begin  with  the  general  equation  for  an  interacting  Hamiltonian. 


H=  ¥(r)Hli>(r)dr  +  -  /  - V)^(f)^(P) 


(Dl) 


Then  we  make  the  following  substitutions  into  our  general  interacting  Hamiltonian. 


e^dr 


(D2) 


*,(0  =  E^(2)^s 


e~ik?d\ 

k 


(D3) 


1  -2 

X{z)  =  X*{z)  =  — — e^5 
VC7T4 


(D4) 


V(r-r' )  =  J2v(p) 


LD 


3—  ip-(r— rf) 


(05) 


We  begin  with  calculating  the  simpler  non-interacting  component  and  then  follow  with 
the  interacting  portion. 


=  [  dfC(z)ie-i5p‘(-lv*  -  U(z))^x(^q-  (D6) 

Hi  =  E  p  /  drx^e-^e  ~~  +  (D7) 

kq 

The  Hamiltonian  operator  within  the  non-interacting  portion  provides  an  average  kinetic 
energy  summed  over  all  particles.  This  gives  us  the  total  kinetic  energy  of  the  system.  The 
potential  portion  is  the  same  as  for  a  harmonic  oscillator  and  since  xiz)  is  the  ground  state 
of  the  system,  the  potential  is  simply  \hio.  Since  we  set  the  condition  that  h  =  m  =  u  =  1 
as  an  assumption  to  simplify  calculations,  this  leaves  a  potential  of  | 
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1  v-  fk2  1 


Hi  =  i-2  £  (  V  +  i  )  /  x(z)e-^-^dp  I  x2(z)dz  did, 

kq 


L2^\  2  '2  II  /  A*  v~'  fc  9 


(D8) 


*  =  £  V  + 


fc2  1 


't: 


i  -  araff 
2  2/  k  q 


(D9) 


Since  \  is  a  constant,  it  can  be  ignored  since  it  affects  every  state  equally.  Therefore,  our 
final  non  interacting  component  is  the  exact  same  as  our  3-D  case. 


2  akaQ 


(DIO) 


Now  we  continue  on  to  the  more  interesting  2-D  contact  interaction  term. 


Hint  =  -  &(i^&(P)V{r-r')V(r)V{r') 


(Dll) 


1  ‘  f  drd?x2(z)x2(z)^  E  V  (p)  eik'pel^'pe^ik'  'pe^H  'pd^d^gdl-;  (D12) 

p 


1  f  f  drd?x2{z)x2{zl)^^v{P)e-i^-P'\ip^-g-^^e-ip^z-z')d^\;a\  (D13) 


P2D 


Pz 


dzdz'nz(z)nz(z)—5(k  +  q  —  p)5(p'  —  k'  —  q') 


dk 

~2n 


—  V{p)e  ipz{z  ^d^dfdldl 


(D14) 


Hint  =  -^j-  I  dze  iPznz(z )  f  dz'ew*'nz(z')  f  T^V(p)e  iPz(z  dp+^^iat  (D15) 


Here  we  use  the  Fourier  transform  (or  momentum  space  equivalent)  nz(kz )  = 

f  dzetkz'zx2(z )  in  order  to  change  the  equation  into  2D  momentum  space. 


Hint  —  2£  f  2^  ^  (p)nz(kz)nz(kz)dp^^^_^dldl 


(D16) 


Hint  2  ®P+rHk-q®]^pV<l2  D  (?) 


(D17) 


kqp 
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Appendix  E:  Derivation  of  Quasi-2D  Dipolar  Interaction  Potential 

The  difference  between  the  contact  and  dipolar  interaction  potential  lies  only  in  the  the 
potential  term.  Dipolar  interactions  follow  the  interaction  potential  Vd(x)  =  d 2 1~3^°s2  - .  To 
approach  this  problem,  we  will  use  spherical  harmonics  and  Bessel  functions. 


H,nt  —  -  /  T^(r)T^(r')D(r  —  r')4>(r)4>(r')drdr' 


(El) 


Below  is  the  potential  for  dipolar  interactions 


Vd(x)  =  d 


,2 1  —  3  cos 2  (9) 
x3 


We  expand  the  potential  using  spherical  harmonics. 


(E2) 


U(f)  =  Aj4yfy20  (E3) 

We  then  Fourier  transform  and  spherically  expand  our  result  into  a  3D  momentum  space 
potential.  What  remains  to  be  done  is  to  convert  our  3D  potential  into  a  quasi  2D  potential. 


=  4tt  £  ilYl^n(k)ji(kx)Yim(x) 


(E4) 


l,m 


Vd(k)  =  7rd2t(3cos2(0fc)  -  1)  where  cos(0k)  =  F 


(E5) 


~  f  dk  ~  ^ 

V,2Dd(k)  =  /  -^n\(k,)Vd(k) 


(E6) 


After  substituting  our  value  for  14 (k)  into  the  equation  we  calculated  above  in  Appendix 
A,  we  use  Mathematica  to  calculate  Vq2Dd(k). 


Vq2Dd{k )  =  fl  (2  -  3a Jnqeq\r 
a/27 tL 


h  1 

repU~ 

where  q  =  -*-=■ 

72 


(E7) 
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(jd  is  a  substitution  used  to  give  the  equation  a  similar  form  regardless  of  our  interaction 
type.  By  making  =  lnhJldd ,  14 (k)  inhabits  the  same  form  whether  in  the  2D  contact  case 
or  the  2D  dipolar  case.  In  this  equation,  acid  =  ^r.  Vq2Dd{k)  has  a  crucial  distinction  between 
it  and  the  2D  contact  case.  The  momentum  dependence  coinciding  with  dipolar  interactions 
makes  calculations  significantly  more  difficult,  than  with  momentum  independent  potentials. 
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Appendix  F:  Solving  the  Hartree  Fock  Bogoliubov  Equations 

Through  our  mathematical  framework,  we  find  a  collection  of  equations  we  call  the 
Hartree  Fock  Bogoliubov  (HFB)  equations. 


h 

£11  (k) 

£20(£) 

“ k 

uk 

vk 

n(k) 

k 


nV(°)  +  Jp^2v(^)n{k) 


nV (0)  +  n0V (k)  +  ^  n(q)V (q  —  k) 


o 


n0V(k ) 


e%- n  +  Y,u(k))  -  (E20(A:)) 


\ 


1  f  ek  ~  V  +  En(fc) 


+  1 


Wi 


1  /  —  F  +  2d ii  (A;) 


CUr 


A 


(ui  +  uI)  k  +  v, 

1 

e^s  -  1 


-  1 


n0 


dk 

pF 


n(fc) 


(Fl) 

(F2) 

(F3) 

(F4) 

(F5) 

(F6) 

(F7) 

(F8) 

(F9) 


We  have  nine  equations  and  nine  unknowns!  In  order  to  calculate  a  self  consistent  solu¬ 
tion,  we  must  ensure  that  particle  number  remains  constant.  We  first  guess  a  solution  for 
the  condensate  fraction  no-  We  use  this  value  to  calculate  the  self  energies  En  and  S2o-  We 
use  the  self  energies  to  calculate  c Uj:,  and  v%.  Once  we  have  the  value  for  we  plug 
it  into  the  Bose  Einstein  distribution  to  find  /g  which  we  use  to  End  our  thermal  density, 
n(k).  We  use  the  thermal  density  to  hnd  the  number  of  excited  atoms  in  the  system. 

However,  we  often  find  the  number  of  atoms  in  our  system,  excited  atoms  plus  condensate 
atoms,  exceed  the  total  particle  number.  In  order  to  fix  this,  we  correct  our  initial  guess  of 
the  condensate  fraction  and  iterate  the  process  using  the  Matlab  code  shown  in  Appendix 
G.  Through  several  iterations,  we  find  a  converged  solution.  We  perform  this  analysis  to 
calculate  all  of  our  results. 
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Appendix  G:  Matlab  Code  for  Solving  Hartree  Fock  Bogoliubov  Equations 

Here  we  include  the  Matlab  code  used  for  solving  our  HFB  equations, 
function  data  =  HFB_Popov_2D_Trident (g,gd,T,ntotal , nOguess , nTguess ,kmin,kmax,klength, ite 

%  HFB  solver 

kvec  =  linspace(kmin,kmax,klength) ; 

nO  =  nOguess; 

if  isempty (nTguess) 
nT  =  zeros (size (kvec) ) ; 
else 

nT  =  nTguess; 
end 

nTold  =  nT; 

maxiter  =  200; 
exitflag  =  0; 
dn0=l ; 
iter  =  1; 

tic 

while  (dn0>le-16  &&  ~exitflag  &&  iter<maxiter) 

if  iter<itertol 
nTfrac  =  iter/itertol ; 
else 

nTfrac  =  1; 


end 
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[wk , uk , vk , SI 1 , S20 , S20T0 , nTnew , mu]  =  BdG (kvec , g , gd , nO , ntotal , nT , T) ; 

if  isreal(wk) 

nT  =  nTfrac*nTnew+(l-nTfrac)*nTold; 

nTtot  =  2*pi/(2*pi)~2*trapz(kvec,kvec.*nT) ; 

nOold  =  nO; 

nO  =  ntotal-nTtot ; 

nO-nOold; 

dnO  =  abs (nO-nOold) ; 

else 

n0=0 ; 

exitflag  =  1; 
end 

nTold  =  nT; 
iter  =  iter+1; 
t  =  toe; 

end 

°/„  find  superfluid  fraction 
nSF  =  SFfun(T, kvec, wk, ntotal) ; 

°/„  find  local  minimum  in  wk 
dwk  =  gradient (wk, kvec) ; 

wkind  =  min(find(gradient(dwk./abs(dwk))==l)) ; 

if  isempty (wkind) 

kroton  =  nan; 

wkroton  =  nan; 

else 

kroton  =  kvec (wkind); 
wkroton  =  wk (wkind) ; 
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end 

fprintf  ( ’nO  =  °/0G  nSF  =  %G  mu  =  "/KAn'  ,nO,nSF,mu) 

data.kroton  =  kroton; 
data.wkroton  =  wkroton; 
data.nSF  =  nSF ; 
data. mu  =  mu; 
data.nO  =  nO; 
data.nTtot  =  nTtot; 
data.nT  =  nT; 
data.n  =  ntotal; 
data.kvec  =  kvec; 
data.wk  =  wk; 
data.uk  =  uk; 
data.vk  =  vk; 
data. Sll  =  Sll; 
data.S20  =  S20; 
data. S20T0  =  S20T0; 
data. iter  =  iter-1; 
data. exitf lag  =  exitflag; 
data.t  =  t; 

%  end  of  HFB.m 

function  [wk,uk,vk,Sll,S20,S20T0,nT,mu]  =  BdG (kvec, g,gd,nO, ntotal, nT,T) 
[Sll ,  S20,S20T0,mu]  =  SEf un (kvec , g, gd,nO , ntotal ,nT) ; 

Ak  =  kvec . ~2/2-mu+Sll ; 

Bk  =  S20 ; 


wk  =  sqrt (Ak . ~2-Bk . ~2) ; 

uk  =  1 . /sqrt ( 1— ( (wk-Ak) . /Bk) . “2) ; 
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vk  =  ( (wk-Ak) . /Bk) . /sqrt (1- ( (wk-Ak) . /Bk) . ~2) ; 

if  T>le-6 

fk  =  1 . /(exp(wk/T)-l) ; 

else 

fk  =  0; 

end 

nT  =  1/2* (uk+vk) .  ~2 . *f k+vk . ~2 ; 

return 

function  [Sll,S20,S20T0,mu]  =  SEfun(kvec,g,gd,nO,ntotal,nT) 
thetavec  =  linspace(0,pi , 128) ; 

Vhn  =  (2*pi/(2*pi)~2)*trapz(kvec,kvec.*Vkfun(kvec,g,gd) . *nT) ; 

Vfn  =  zeros (size (kvec) ) ; 
for  ii  =  1 : length (kvec) 
k  =  kvec(ii) ; 

Vfvec  =  zeros (size (kvec) ) ; 
for  jj  =  1 : length (kvec) 
q  =  kvec(j j) ; 

Vfvec(j j)  =  2*trapz (thetavec, Vkfun (sqrt (q. ~2+k. ~2-2*k*q*cos (thetavec) ) ,g,gd))  ; 
end 

Vfn(ii)  =  1/ (2*pi) ~2*trapz (kvec, kvec . *Vfvec . *nT) ; 
end 

Sll  =  ntotal*Vkfun(0,g,gd)+nO*Vkfun(kvec,g,gd)+Vfn; 

S20  =  nO*Vkfun(kvec,g,gd) ; 

S20T0  =  nO*Vkfun(kvec,g,gd) ; 
mu  =  ntotal*Vkfun(0,g,gd)+Vhn; 
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return 

function  Vk  =  Vkfun(kin,g,gd) 

k  =  kin( : )/sqrt(2) ; 
ks  =  k(k<=25); 
kl  =  k(k>25) ; 

vals  =  exp(ks . ~2) . *erf c (ks) ; 

vail  =  0 . 0763957775554068*exp (-0 . 049875774108393764*kl) ; 

Vks  =  g+gd*(2-3*sqrt(pi)*ks . *vals) ; 

Vkl  =  g+gd* (vail . /kl-1) ; 

Vk  =  [Vks ; Vkl] ; 

Vk  =  reshape (Vk, size (kin) ) ; 

return 

function  nSF  =  SFfun(T,kvec,wk,ntotal) 

nSF  =  ntotal+l/2*2*pi/ (2*pi) ~2*trapz(kvec ,kvec . ~3 . * (1 . / (l-cosh(wk/T) ) )/2/T) ; 
if  nSF<0 
nSF=0 ; 
end 


return 


